
    ##########################################################################################
    ### R code for Figure 4 of Supplemental Material
    ### : Monte Carlo Experiment to Show that Interest Group Ratings tend to overemphasize polarization:
    ### : Multilple Iterations
    ##########################################################################################
    
    library(mvtnorm)
    library(MCMCpack)

    ### The First Scenario
    ##########################################################################################
    
    ### I. Data Generating Process:
    ### 1. Fix the numbers of legislators and roll call votes.

    M <- 100           # Number of roll call votes
    n.D <- 258
    n.R <- 176
    N <- n.D + n.R     # Number of Legislators

    ### Number of Simulation
    G <- 50

    Ideal <- IRT   <- matrix(NA, N, G)
    ADA1 <- ADA2  <- matrix(NA, N, G)
    n.close.votes1 <- rep(NA, G)
    n.close.votes2 <- rep(NA, G)

    Cutoff <- matrix(NA, M, G)
    n.close.votes1 <- n.close.votes2 <- rep(NA, G)
    
    ### LOOP over simulations #############
    for (g in 1:G){

    # Generate Ideal Points
    Ideal[,g] <- c(rnorm(n.D, mean=-1, sd=1), rnorm(n.R, mean=1, sd=1))
    Ideal[c(1,N),g] <- c(-2.5, 2.5) # Ideal points for constraints

    # Generate Cut points
    Cutoff[,g] <- rnorm(M, mean=0, sd=1)
    
    # Generate voting data
    Y <- matrix(NA, N, M)
    row.names(Y) <- paste("L", 1:N, sep="")
    # the following code is same as "Y ~ dbern(P)" in WinBUGS
    for (i in 1:N){
    for (j in 1:M){
    if (Ideal[i,g] >= Cutoff[j,g]) Y[i,j] <- 0.0
    else Y[i,j] <- 1.0
    }}

    # Compute vote margins
    Yeas <- matrix(NA, M, 1)
    for(j in 1:M){
    Yeas[j,1] <- sum(Y[,j])/N
    }

    # Select votes with close margins (40-60%)
    close.votes1 <- rep(NA, M)
    for (j in 1:M){    
    if( Yeas[j,1] > .4 & Yeas[j,1] <= .6) close.votes1[j] <- 1.0
    else close.votes1[j] <- 0.0
    }

    Y.close1 <- Y[, close.votes1==1]
    n.close.votes1[g] <- n.Y.close1 <- ncol(Y.close1)

    # Select votes with close margins (20-80%)
    close.votes2 <- rep(NA, M)
    for (j in 1:M){    
    if( Yeas[j,1] > .2 & Yeas[j,1] <= .8) close.votes2[j] <- 1.0
    else close.votes2[j] <- 0.0
    }

    Y.close2 <- Y[, close.votes2==1]
    n.close.votes2[g] <- n.Y.close2 <- ncol(Y.close2)
  
    ## Estimation
    # First, IRT estimation using MCMCirt1d
                                         
    MCMC <- MCMCirt1d(Y, theta.constraints=list(L1="-", L434="+"),  
                      burnin=1000, mcmc=2000, thin=2, store.item=FALSE, store.ability=TRUE)   
    iter <- nrow(MCMC)

    # Save the mean ideal point estimates
    IRT[,g] <- apply(MCMC[,1:N], 2, FUN=mean)
  
    ## Second, ADA score with votes (40-60% margins)
    ADA.row1  <- apply(Y.close1, 1, FUN=sum)
    ADA1[,g] <- 100*(ADA.row1/n.Y.close1)

    ## Third, ADA score with votes (20-80% margins)
    ADA.row2  <- apply(Y.close2, 1, FUN=sum)
    ADA2[,g] <- 100*(ADA.row2/n.Y.close2) 
    }

    ### End of Loop ##########################

    # Rewrite the names to avoid being overwritten
    MC1_Ideal <- Ideal 
    MC1_IRT   <- IRT
    MC1_ADA1  <- ADA1
    MC1_ADA2  <- ADA2 



    ### The Second Scenario
    ##########################################################################################

    ### I. Data Generating Process:
    # 1. Fix the numbers of legislators and roll call votes.

    M <- 100           # Number of roll call votes
    n.D <- 258
    n.R <- 176
    N <- n.D + n.R     # Number of Legislators

    Ideal <- IRT   <- matrix(NA, N, G)
    ADA1 <- ADA2  <- matrix(NA, N, G)
    n.close.votes1 <- rep(NA, G)
    n.close.votes2 <- rep(NA, G)

    Cutoff <- matrix(NA, M, G)
    n.close.votes1 <- n.close.votes2 <- rep(NA, G)
    
    
    ### LOOP over simulations ########################
    for (g in 1:G){
    # Generate Ideal Points
    Ideal[,g] <- c(rnorm(n.D, mean=-2, sd=1), rnorm(n.R, mean=2, sd=1))
    Ideal[c(1,N),g] <- c(-2.5, 2.5) # Ideal points for constraints

    # Generate Cut points
    Cutoff[,g] <- rnorm(M, mean=0, sd=1)
    
    # Generate voting data
    Y <- matrix(NA, N, M)
    row.names(Y) <- paste("L", 1:N, sep="")
    # the following code is same as "Y ~ dbern(P)" in WinBUGS
    for (i in 1:N){
    for (j in 1:M){
    if (Ideal[i,g] >= Cutoff[j,g]) Y[i,j] <- 0.0
    else Y[i,j] <- 1.0
    }}

    # Compute vote margins
    Yeas <- matrix(NA, M, 1)
    for(j in 1:M){
    Yeas[j,1] <- sum(Y[,j])/N
    }

    # Select votes with close margins (40-60%)
    close.votes1 <- rep(NA, M)
    for (j in 1:M){    
    if( Yeas[j,1] > .4 & Yeas[j,1] <= .6) close.votes1[j] <- 1.0
    else close.votes1[j] <- 0.0
    }

    Y.close1 <- Y[, close.votes1==1]
    n.close.votes1[g] <- n.Y.close1 <- ncol(Y.close1)

    # Select votes with close margins (20-80%)
    close.votes2 <- rep(NA, M)
    for (j in 1:M){    
    if( Yeas[j,1] > .2 & Yeas[j,1] <= .8) close.votes2[j] <- 1.0
    else close.votes2[j] <- 0.0
    }

    Y.close2 <- Y[, close.votes2==1]
    n.close.votes2[g] <- n.Y.close2 <- ncol(Y.close2)
 
    ## Estimation
    # First, IRT estimation using MCMCirt1d
                                         
    MCMC <- MCMCirt1d(Y, theta.constraints=list(L1="-", L434="+"),  
                      burnin=1000, mcmc=2000, thin=2, store.item=FALSE, store.ability=TRUE)   
    iter <- nrow(MCMC)

    # Save the mean ideal point estimates
    IRT[,g] <- apply(MCMC[,1:N], 2, FUN=mean)
  
    ## Second, ADA score with votes (40-60% margins)
    ADA.row1  <- apply(Y.close1, 1, FUN=sum)
    ADA1[,g] <- 100*(ADA.row1/n.Y.close1)

    ## Third, ADA score with votes (20-80% margins)
    ADA.row2  <- apply(Y.close2, 1, FUN=sum)
    ADA2[,g] <- 100*(ADA.row2/n.Y.close2) 
    }

    ### End of Loop ##############################################

    # Rewrite the names to avoid being overwritten
    MC2_Ideal <- Ideal 
    MC2_IRT   <- IRT
    MC2_ADA1  <- ADA1
    MC2_ADA2  <- ADA2 

    # Reverse the direct of ADA scores to make them consistent with IRT

    MC1_ADA1 <- - MC1_ADA1
    MC1_ADA2 <- - MC1_ADA2
    
    MC2_ADA1 <- - MC2_ADA1
    MC2_ADA2 <- - MC2_ADA2

    postscript("Supplemental_Figure4.eps", height=11, width=6.5, horizontal=F)
    par(mfcol=c(3, 2))
    plot(MC1_Ideal[,1], MC1_IRT[,1], type='p', col='grey', 
         xlim=c(min(MC1_Ideal), max(MC1_Ideal)), 
         ylim=c(min(MC1_IRT), max(MC1_IRT)), 
         xlab="True Ideal Points", ylab="IRT Estimate")
    title(main="DEM ~ N(-1,1); GOP ~ N(1,1)", font.main=1)
    #segments(min(MC1_Ideal[,1]), min(MC1_IRT[,1]), max(MC1_Ideal[,1]), max(MC1_IRT[,1]), lty=1)
    for(g in 1:G){
    points(MC1_Ideal[,g], MC1_IRT[,g], type='p', col='grey') 
    }
    plot(MC1_Ideal[,1], MC1_ADA1[,1], type='p', col='grey', 
         xlim=c(min(MC1_Ideal), max(MC1_Ideal)), 
         ylim=c(min(MC1_ADA1[,1]), max(MC1_ADA1[,1])), 
         xlab="True Ideal Points", ylab="ADA (Majority Size = 60% or less)")
    title(main="DEM ~ N(-1,1); GOP ~ N(1,1)", font.main=1)
    #segments(min(MC1_Ideal[,g]), min(ADA[,g]), max(MC1_Ideal[,g]), max(ADA[,g]), lty=1)
    for(g in 1:G){
    points(MC1_Ideal[,g], MC1_ADA1[,g], type='p', col='grey') 
    }
    plot(MC1_Ideal[,1], MC1_ADA2[,1], type='p', col='grey', 
         xlim=c(min(MC1_Ideal), max(MC1_Ideal)), 
         ylim=c(min(MC1_ADA2[,1]), max(MC1_ADA2[,1])), 
         xlab="True Ideal Points", ylab="ADA (Majority Size = 80% or less)")
    #segments(min(MC1_Ideal[,g]), min(MC1_ADA2[,g]), max(MC1_Ideal[,g]), max(MC1_ADA2[,g]), lty=1)
    title(main="DEM ~ N(-1,1); GOP ~ N(1,1)", font.main=1)
    for(g in 1:G){
    points(MC1_Ideal[,g], MC1_ADA2[,g], type='p', col='grey') 
    }

    plot(MC2_Ideal[,1], MC2_IRT[,1], type='p', col='grey', 
         xlim=c(min(MC2_Ideal), max(MC2_Ideal)), 
         ylim=c(min(MC2_IRT), max(MC2_IRT)), 
         xlab="True Ideal Points", ylab="IRT Estimate")
    #segments(min(MC2_Ideal[,1]), min(MC2_IRT[,1]), max(MC2_Ideal[,1]), max(MC2_IRT[,1]), lty=1)
    title(main="DEM ~ N(-2,1); GOP ~ N(2,1)", font.main=1)
    for(g in 1:G){
    points(MC2_Ideal[,g], MC2_IRT[,g], type='p', col='grey') 
    }
    plot(MC2_Ideal[,1], MC2_ADA1[,1], type='p', col='grey', 
         xlim=c(min(MC2_Ideal), max(MC2_Ideal)), 
         ylim=c(min(MC2_ADA1[,1]), max(MC2_ADA1[,1])), 
         xlab="True Ideal Points", ylab="ADA (Majority Size = 60% or less)")
    #segments(min(MC2_Ideal[,g]), min(ADA[,g]), max(MC2_Ideal[,g]), max(ADA[,g]), lty=1)
    title(main="DEM ~ N(-2,1); GOP ~ N(2,1)", font.main=1)
    for(g in 1:G){
    points(MC2_Ideal[,g], MC2_ADA1[,g], type='p', col='grey') 
    }
    plot(MC2_Ideal[,1], MC2_ADA2[,1], type='p', col='grey', 
         xlim=c(min(MC2_Ideal), max(MC2_Ideal)), 
         ylim=c(min(MC2_ADA2[,1]), max(MC2_ADA2[,1])), 
         xlab="True Ideal Points", ylab="ADA (Majority Size = 80% or less)")
    #segments(min(MC2_Ideal[,g]), min(MC2_ADA2[,g]), max(MC2_Ideal[,g]), max(MC2_ADA2[,g]), lty=1)
    title(main="DEM ~ N(-2,1); GOP ~ N(2,1)", font.main=1)
    for(g in 1:G){
    points(MC2_Ideal[,g], MC2_ADA2[,g], type='p', col='grey') 
    }

    dev.off()    
    
    
    
